#include "header.h"

/* -------------------------------------------------------------------------- */
void xyzflux(double *var,double *Xflux,double* Yflux, double* Zflux){
/* -------------------------------------------------------------------------- */
double vel2,ener;

vel2	= var[1]*var[1] + var[2]*var[2] + var[3]*var[3];
ener    = var[4]/(GAMMA-1.0) + 0.5*var[0]*vel2;

Xflux[0] = var[0]*var[1];
Xflux[1] = var[4] + var[0]*var[1]*var[1];
Xflux[2] = var[0]*var[1]*var[2];
Xflux[3] = var[0]*var[1]*var[3];
Xflux[4] = (ener + var[4])*var[1];

Yflux[0] = var[0]*var[2];
Yflux[1] = var[0]*var[1]*var[2];
Yflux[2] = var[4] + var[0]*var[2]*var[2];
Yflux[3] = var[0]*var[2]*var[3];
Yflux[4] = (ener + var[4])*var[2];

Zflux[0] = var[0]*var[3];
Zflux[1] = var[0]*var[1]*var[3];
Zflux[2] = var[0]*var[2]*var[3];
Zflux[3] = var[4] + var[0]*var[3]*var[3];
Zflux[4] = (ener + var[4])*var[3];
}
/* -------------------------------------------------------------------------- */
void FluxKinetic(double* vleft, double* vright, double* flux){
/* -------------------------------------------------------------------------- */
int i;
double fp[5],fm[5];
xflux_plus(vleft,fp);
xflux_minus(vright,fm);
for(i=0;i<5;i++) flux[i] = fp[i] + fm[i];
}
/* -------------------------------------------------------------------------- */
void xflux_plus(double *var,double *flux){
/* -------------------------------------------------------------------------- */
double beta,mach,vel2,err_fun,exp_fun,fact,ener;

beta 	= 0.5*var[0]/var[4];
mach	= var[1]*sqrt(beta);
err_fun	= 0.5*(1.0 + ERRF(mach));
exp_fun	= 0.5*exp(-mach*mach)/sqrt(M_PI*beta);
fact	= var[1]*err_fun + exp_fun;
vel2	= var[1]*var[1] + var[2]*var[2] + var[3]*var[3];
ener    = var[4]/(GAMMA-1.0) + 0.5*var[0]*vel2;

flux[0]	= var[0]*fact;
flux[1]	= (var[4]+var[0]*var[1]*var[1])*err_fun + var[0]*var[1]*exp_fun;
flux[2] = var[0]*var[2]*fact;
flux[3] = var[0]*var[3]*fact;
flux[4] = var[1]*(ener + var[4])*err_fun + (ener + 0.5*var[4])*exp_fun;
}
/* -------------------------------------------------------------------------- */
void xflux_minus(double *var,double *flux){
/* -------------------------------------------------------------------------- */
double beta,mach,vel2,err_fun,exp_fun,fact,ener;

beta 	= 0.5*var[0]/var[4];
mach	= var[1]*sqrt(beta);
err_fun	= 0.5*(1.0 - ERRF(mach));
exp_fun	= 0.5*exp(-mach*mach)/sqrt(M_PI*beta);
fact	= var[1]*err_fun - exp_fun;
vel2	= var[1]*var[1] + var[2]*var[2] + var[3]*var[3];
ener    = var[4]/(GAMMA-1.0) + 0.5*var[0]*vel2;

flux[0]	= var[0]*fact;
flux[1]	= (var[4]+var[0]*var[1]*var[1])*err_fun - var[0]*var[1]*exp_fun;
flux[2] = var[0]*var[2]*fact;
flux[3] = var[0]*var[3]*fact;
flux[4] = var[1]*(ener + var[4])*err_fun - (ener + 0.5*var[4])*exp_fun;
}
/* -------------------------------------------------------------------------- */
void yflux_plus(double *var,double *flux){
/* -------------------------------------------------------------------------- */
double beta,mach,vel2,err_fun,exp_fun,fact,ener;

beta    = 0.5*var[0]/var[4];
mach    = var[2]*sqrt(beta);
err_fun = 0.5*(1.0 + erf(mach));
exp_fun = 0.5*exp(-mach*mach)/sqrt(M_PI*beta);
fact    = var[2]*err_fun + exp_fun;
vel2    = var[1]*var[1] + var[2]*var[2] + var[3]*var[3];
ener    = var[4]/(GAMMA-1.0) + 0.5*var[0]*vel2;

flux[0] = var[0]*fact;
flux[1] = var[0]*var[1]*fact;
flux[2] = (var[4] + var[0]*var[2]*var[2])*err_fun + var[0]*var[2]*exp_fun;
flux[3] = var[0]*var[3]*fact;
flux[4] = var[2]*(var[4] + ener)*err_fun + (0.5*var[4] + ener)*exp_fun;
}
/* -------------------------------------------------------------------------- */
void yflux_minus(double *var,double *flux){
/* -------------------------------------------------------------------------- */
double beta,mach,vel2,err_fun,exp_fun,fact,ener;
 
beta    = 0.5*var[0]/var[4];
mach    = var[2]*sqrt(beta);
err_fun = 0.5*(1.0 - erf(mach));
exp_fun = 0.5*exp(-mach*mach)/sqrt(M_PI*beta);
fact    = var[2]*err_fun - exp_fun;
vel2    = var[1]*var[1] + var[2]*var[2] + var[3]*var[3]; 
ener    = var[4]/(GAMMA-1.0) + 0.5*var[0]*vel2;
 
flux[0] = var[0]*fact;
flux[1] = var[0]*var[1]*fact;
flux[2] = (var[4] + var[0]*var[2]*var[2])*err_fun - var[0]*var[2]*exp_fun;
flux[3] = var[0]*var[3]*fact;
flux[4] = var[2]*(var[4] + ener)*err_fun - (0.5*var[4] + ener)*exp_fun;
}
/* -------------------------------------------------------------------------- */
void yfluxN_plus(double *var,double *flux){
/* -------------------------------------------------------------------------- */
double beta,mach,vel2,err_fun,exp_fun,fact;
double ener;

beta 	= 0.5*var[0]/var[4];
mach	= var[1]*sqrt(beta);
err_fun	= 0.5*(1.0 + ERRF(mach));
exp_fun	= 0.5*exp(-mach*mach)/sqrt(M_PI*beta);
fact	= var[1]*err_fun + exp_fun;
vel2	= var[1]*var[1] + var[2]*var[2] + var[3]*var[3];
ener    = var[4]/(GAMMA-1.0) + 0.5*var[0]*vel2;

flux[0]	= var[0]*var[2]*err_fun;
flux[1]	= var[0]*var[2]*fact;
flux[2] = (var[4] + var[0]*var[2]*var[2])*err_fun;
flux[3] = var[0]*var[2]*var[3]*err_fun;
flux[4] = var[2]*( var[4] + ener)*err_fun + 0.5*var[0]*var[1]*var[2]*exp_fun;
}
/* -------------------------------------------------------------------------- */
void yfluxN_minus(double *var,double *flux){
/* -------------------------------------------------------------------------- */
double beta,mach,vel2,err_fun,exp_fun,fact;
double ener;

beta 	= 0.5*var[0]/var[4];
mach	= var[1]*sqrt(beta);
err_fun	= 0.5*(1.0 - ERRF(mach));
exp_fun	= 0.5*exp(-mach*mach)/sqrt(M_PI*beta);
fact	= var[1]*err_fun - exp_fun;
vel2	= var[1]*var[1] + var[2]*var[2] + var[3]*var[3];
ener    = var[4]/(GAMMA-1.0) + 0.5*var[0]*vel2;

flux[0]	= var[0]*var[2]*err_fun;
flux[1]	= var[0]*var[2]*fact;
flux[2] = (var[4] + var[0]*var[2]*var[2])*err_fun;
flux[3] = var[0]*var[2]*var[3]*err_fun;
flux[4] = var[2]*( var[4] + ener)*err_fun - 0.5*var[0]*var[1]*var[2]*exp_fun;
}
/* -------------------------------------------------------------------------- */
void zfluxN_plus(double *var,double *flux){
/* -------------------------------------------------------------------------- */
double beta,mach,vel2,err_fun,exp_fun,fact;
double ener;

beta    = 0.5*var[0]/var[4];
mach    = var[1]*sqrt(beta);
err_fun = 0.5*(1.0 + ERRF(mach));
exp_fun = 0.5*exp(-mach*mach)/sqrt(M_PI*beta);
fact    = var[1]*err_fun + exp_fun;
vel2    = var[1]*var[1] + var[2]*var[2] + var[3]*var[3];
ener    = var[4]/(GAMMA-1.0) + 0.5*var[0]*vel2;

flux[0] = var[0]*var[3]*err_fun;
flux[1] = var[0]*var[3]*fact;
flux[2] = var[0]*var[3]*var[2]*err_fun;
flux[3] = (var[4] + var[0]*var[3]*var[3])*err_fun;
flux[4] = var[3]*( var[4] + ener)*err_fun + 0.5*var[0]*var[1]*var[3]*exp_fun;
}
/* -------------------------------------------------------------------------- */
void zfluxN_minus(double *var,double *flux){
/* -------------------------------------------------------------------------- */
double beta,mach,vel2,err_fun,exp_fun,fact;
double ener;

beta    = 0.5*var[0]/var[4];
mach    = var[1]*sqrt(beta);
err_fun = 0.5*(1.0 - ERRF(mach));
exp_fun = 0.5*exp(-mach*mach)/sqrt(M_PI*beta);
fact    = var[1]*err_fun - exp_fun;
vel2    = var[1]*var[1] + var[2]*var[2] + var[3]*var[3];
ener    = var[4]/(GAMMA-1.0) + 0.5*var[0]*vel2;

flux[0] = var[0]*var[3]*err_fun;
flux[1] = var[0]*var[3]*fact;
flux[2] = var[0]*var[3]*var[2]*err_fun;
flux[3] = (var[4] + var[0]*var[3]*var[3])*err_fun;
flux[4] = var[3]*( var[4] + ener)*err_fun - 0.5*var[0]*var[1]*var[3]*exp_fun;
}
/* -------------------------------------------------------------------------- */
